home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / External Modules / External modules sources / Pascal / FFT.p next >
Text File  |  1994-09-28  |  11KB  |  396 lines
  1. { This unit contains the function 'fft' which, together with 350 additional }
  2. { programs useful in scientific work, is found in the book 'Numerical Recipes: }
  3. { The Art of Scientific Computing', available from Cambridge University Press. }
  4. { It is used here by permission of the authors. }
  5.  
  6. {************************************************************************************}
  7. {  FFT.p                                                                      }
  8. {                                                                           }
  9. {                                                                           }
  10. {  Version 26.9.94                                                      }
  11. {************************************************************************************}
  12.  
  13.  
  14. unit user;
  15.  
  16. interface
  17.     uses
  18.         proFit_interface;
  19.  
  20.     procedure main (selector: integer; pb: ExtModulesParamBlockPtr);
  21.  
  22.  
  23.  
  24. implementation
  25.  
  26. { note: MPW users must make sure that the procedure main is at the beginning of the compiled code }
  27. { under Think Pascal, this is cared for by the compiler }
  28. { We let main call a function mainMain to make sure that the code starts with a jump to }
  29. { our entry point even when compiling under MPW Pascal }
  30.  
  31.     procedure mainMain (selector: integer; pb: ExtModulesParamBlockPtr);
  32.     forward;
  33.     procedure main (selector: integer; pb: ExtModulesParamBlockPtr);
  34.  
  35.     begin
  36.         mainMain(selector, pb);
  37.     end;
  38.  
  39.  
  40.  
  41.  
  42. {--------------------------------------------------------------------------------------------------------}
  43.  
  44.  
  45.  
  46.     type
  47.         data = array[1..8192] of extended;
  48.         dataptr = ^data;
  49.  
  50.  
  51.     procedure WriteStringAndNumber (s: str255; num: longint);
  52.     begin
  53.         WriteString(s);
  54.         NumToString(num, s);
  55.         WritelnString(s);
  56.     end;        { WriteStringAndNumber }
  57.  
  58.  
  59.     procedure fft (x: dataptr; num: integer; inverse: boolean);
  60.         var
  61.             ii, jj, n, mmax, m, j, istep, i: integer;
  62.             wtemp, wr, wpr, wpi, wi, theta, tempr, tempi: extended;
  63.     begin
  64.         n := num;
  65.         j := 1;
  66.         for ii := 1 to (n div 2) do
  67.             begin
  68.                 i := 2 * ii - 1;
  69.                 if j > i then
  70.                     begin
  71.                         tempr := x^[j];
  72.                         tempi := x^[j + 1];
  73.                         x^[j] := x^[i];
  74.                         x^[j + 1] := x^[i + 1];
  75.                         x^[i] := tempr;
  76.                         x^[i + 1] := tempi;
  77.                     end;
  78.                 m := n div 2;
  79.                 while ((m >= 2) & (j > m)) do
  80.                     begin
  81.                         j := j - m;
  82.                         m := m div 2;
  83.                     end;
  84.                 j := j + m;
  85.                 if TestStop then        { test if proFit or a user tries to stop the program }
  86.                     exit(fft);
  87.             end;        { for ii := 1 to (n div 2) do }
  88.  
  89.         mmax := 2;
  90.         while n > mmax do
  91.             begin
  92.                 istep := 2 * mmax;
  93.                 theta := 6.283185307959 / mmax;
  94.                 if inverse then
  95.                     theta := -theta;
  96.                 wpr := -2 * sqr(sin(0.5 * theta));
  97.                 wpi := sin(theta);
  98.                 wr := 1;
  99.                 wi := 0;
  100.                 for ii := 1 to (mmax div 2) do
  101.                     begin
  102.                         m := 2 * ii - 1;
  103.                         for jj := 0 to ((n - m) div istep) do
  104.                             begin
  105.                                 i := m + jj * istep;
  106.                                 j := i + mmax;
  107.                                 tempr := wr * x^[j] - wi * x^[j + 1];
  108.                                 tempi := wr * x^[j + 1] + wi * x^[j];
  109.                                 x^[j] := x^[i] - tempr;
  110.                                 x^[j + 1] := x^[i + 1] - tempi;
  111.                                 x^[i] := x^[i] + tempr;
  112.                                 x^[i + 1] := x^[i + 1] + tempi;
  113.                             end;        { for jj := 0 to ((n - m) div istep) do }
  114.                         wtemp := wr;
  115.                         wr := wr + wr * wpr - wi * wpi;
  116.                         wi := wi + wi * wpr + wtemp * wpi;
  117.                         if TestStop then        { test if proFit or a user tries to stop the program }
  118.                             exit(fft);
  119.                     end;        { for ii := 1 to (mmax div 2) do }
  120.                 mmax := istep;
  121.             end;        { while n>mmax }
  122.     end;        { fft }
  123.  
  124. {-----------------------------------------}
  125.  
  126.     procedure shiftFFT (x: dataptr; np: integer);
  127.         var
  128.             xx: extended;
  129.             i, np2: integer;
  130.     begin
  131.         for i := 1 to np do
  132.             begin
  133.                 xx := x^[i];
  134.                 x^[i] := x^[i + np];
  135.                 x^[i + np] := xx;
  136.             end;
  137.     end;    { shiftFFT }
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144. {************************************************************************************}
  145.  
  146.     procedure SetUp (var moduleKind: integer;    { set moduleKind to isFunction or isProgram }
  147.                                     var name: Str255;             { the name of the program or function }
  148.                                     var requiredGlobals: longint;     { the number of bytes to be allocated in ExtModulesParamBlock.globals }
  149.                                         { set requiredGlobals to 0 if you don't use this feature }
  150.                                     pb: ExtModulesParamBlockPtr);    { the complete parameter block passed by pro Fit to the }
  151.                                         { routines defined in this file. In most cases it can be ignored }
  152. { SetUp is called once when the external module is linked to pro Fit }
  153.     begin
  154.         moduleKind := isProgram;
  155.         name := 'Complex FFT';
  156.         requiredGlobals := 0;
  157.     end;
  158.  
  159.  
  160.  
  161.  
  162. {************************************************************************************}
  163.  
  164.     procedure InitializeProg (pb: ExtModulesParamBlockPtr);
  165. { Can be left emtpy if not needed. }
  166. { called when the external module is linked to proFit after SetUp was called }
  167. { can be used to inititialize global variables, etc. }
  168.     begin
  169.         pb^.V[1] := 0; {used as a flag to find out when the  FFT program is called for the first time}
  170.         pb^.V[2] := 64;{ the default input parameters}
  171.         pb^.V[3] := 1;
  172.         pb^.V[4] := 4;
  173.         pb^.V[5] := 0;
  174.         pb^.V[6] := 1;
  175.  
  176.     end;
  177.  
  178.  
  179.  
  180.  
  181.  
  182. {************************************************************************************}
  183.  
  184.     procedure Run (pb: ExtModulesParamBlockPtr);
  185. { pro Fit calls this function when the name of the program is chosen from the }
  186. { Run Program submenu in the menu Calc }
  187.         var
  188.             s1, s2, s3, s4, s5, st: str255;
  189.             r: inputRec;
  190.             ex1, ex2, ex3, ex4, ex5, x0, xinc, y0, yinc, yoff: extended;
  191.             i, xCol, yCol: integer;
  192.             x: dataptr;        { pointer to the data array }
  193.             nrData: integer;    { the number of data (integer power of 2) }
  194.             inverse: boolean;    { false if normal FFT should be done, true for reverse FFT }
  195.             answer: integer;
  196.  
  197.         procedure CleanUpRun;
  198.         begin
  199.             StopExecution;
  200.             if x <> nil then
  201.                 DisposPtr(pointer(x));    { we don't need the memory any more }
  202.             exit(Run);                { terminate this call of run }
  203.         end;
  204.  
  205.     begin
  206.         if pb^.V[1] = 0 then
  207.             begin        {  for the first time this program runs }
  208.                 pb^.V[1] := 1;
  209.                 if AskBox(answer, 'Do you want to print some information about this program to the Results window ?') then
  210.                     exit(Run)                    { stop button pressed }
  211.                 else if answer = 1 then        { yes button pressed }
  212.                     begin
  213.                         WritelnString('The number of data points must be between 4 and  4096.');
  214.                         WritelnString('(The program sets the number of data points to the next smaller power of 2)');
  215.                         WritelnString('"First column No of data" means the column number where the "time" values are stored.');
  216.                         WritelnString('The real and imaginary values of your data must be in the succeeding two columns.');
  217.                         WritelnString('"First column No of FFT" means the column number where the "frequency" values will be stored.');
  218.                         WritelnString('The real and imaginary part of the Fourier transformed data will be in the succeeding two columns.');
  219.                         WritelnString('"Offset of FFT" defines an additive constant for the frequency values.');
  220.                         exit(Run);
  221.                     end;
  222.             end;
  223.  
  224.         ex1 := pb^.v[2];                        { do dialog about dataNr, ColumnNrs }
  225.         s1 := 'Number of data points';
  226.         r[1].x := pointer(@ex1);
  227.         r[1].s := pointer(@s1);
  228.         ex2 := pb^.v[3];
  229.         s2 := 'First column No. of data';
  230.         r[2].x := pointer(@ex2);
  231.         r[2].s := pointer(@s2);
  232.         ex3 := pb^.v[4];
  233.         s3 := 'First column No of FFT';
  234.         r[3].x := pointer(@ex3);
  235.         r[3].s := pointer(@s3);
  236.         ex4 := pb^.v[5];
  237.         s4 := '    Offset of FFT';
  238.         r[4].x := pointer(@ex4);
  239.         r[4].s := pointer(@s4);
  240.         ex5 := pb^.v[6];
  241.         s5 := 'FFT (>0) ; inverse FFT (<0) ?';
  242.         r[5].x := pointer(@ex5);
  243.         r[5].s := pointer(@s5);
  244.         if InputBox(5, r) then
  245.             exit(Run)    { if something is not ok we quit }
  246.         else
  247.             begin    { if everything is ok we do something }
  248.                 ex1 := abs(ex1);
  249.                 if (ex1 > 4096) | (ex1 < 4) then
  250.                     exit(Run);
  251.                 nrData := round(exp(ln(2) * trunc(ln(ex1) / ln(2))));
  252.                 ex2 := abs(ex2);
  253.                 if (ex2 > 97) | (ex2 < 1) then
  254.                     exit(Run);
  255.                 xCol := round(ex2);
  256.                 ex3 := abs(ex3);
  257.                 if (ex3 > 97) | (ex3 < 1) then
  258.                     exit(Run);
  259.                 yCol := round(ex3);
  260.                 yoff := ex4;
  261.                 if ex5 >= 0 then
  262.                     inverse := false
  263.                 else
  264.                     inverse := true;
  265.                 pb^.v[2] := nrData;{set the default parameter to the chosen values, for the next time FFT is called}
  266.                 pb^.v[3] := ex2;
  267.                 pb^.v[4] := ex3;
  268.                 pb^.v[5] := ex4;
  269.                 if inverse then
  270.                     pb^.v[6] := -1
  271.                 else
  272.                     pb^.v[6] := 1;
  273.  
  274.  
  275.  
  276.                 x := dataptr(newPtr(sizeof(extended) * longint(2) * nrData));            { prepare memory for data }
  277.                 if x = nil then
  278.                     exit(run);
  279.  
  280.  
  281.                 if TestData(1, xCol) then
  282.                     begin
  283.                         x0 := GetData(1, xCol);
  284.                         if TestData(2, xCol) then
  285.                             xinc := GetData(2, xCol) - x0;
  286.                     end;
  287.  
  288.                 for i := 1 to NrData do                { read data from window }
  289.                     begin
  290.                         if TestData(i, xCol + 1) then        { test if data value is ok }
  291.                             x^[2 * i - 1] := GetData(i, xCol + 1)
  292.                         else
  293.                             x^[2 * i - 1] := 0;
  294.                         if TestData(i, xCol + 2) then
  295.                             x^[2 * i] := GetData(i, xCol + 2)
  296.                         else
  297.                             x^[2 * i] := 0;
  298.                         if TestStop then
  299.                             CleanUpRun;
  300.                     end;
  301.  
  302.                 fft(x, 2 * NrData, inverse);            { fft calculations }
  303.                 yinc := 1 / NrData / xinc;
  304.                 if not inverse then
  305.                     begin
  306.                         shiftFFT(x, NrData);
  307.                         y0 := yoff - 1 / xinc / 2;
  308.                     end
  309.                 else
  310.                     begin
  311.                         y0 := yoff;
  312.                         for i := 1 to NrData div 2 do
  313.                             begin
  314.                                 x^[4 * i - 3] := x^[4 * i - 3] / NrData;
  315.                                 x^[4 * i - 2] := x^[4 * i - 2] / NrData;
  316.                                 x^[4 * i - 1] := -x^[4 * i - 1] / NrData;
  317.                                 x^[4 * i] := -x^[4 * i] / NrData;
  318.                             end;
  319.                     end;
  320.                 if TestStop then
  321.                     CleanUpRun;
  322.  
  323.                 SetColName(yCol, 'frequency');
  324.                 SetColName(yCol + 1, 'FFT real');
  325.                 SetColName(yCol + 2, 'FFT imag.');
  326.                 for i := 1 to NrData do                    { output of data into window }
  327.                     begin
  328.                         SetData(i, yCol, y0 + (i - 1) * yinc);
  329.                         SetData(i, yCol + 1, x^[2 * i - 1]);
  330.                         SetData(i, yCol + 2, x^[2 * i]);
  331.                         if TestStop then
  332.                             CleanUpRun;
  333.                     end;
  334.  
  335.                 if x <> nil then
  336.                     DisposPtr(Ptr(x));{ we don't need the memory any more }
  337.             end;
  338.     end;{Run}
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346. {************************************************************************************}
  347.  
  348.     procedure CleanUp (pb: ExtModulesParamBlockPtr);
  349.     { called when the function or program is removed from pro Fit's menus }
  350.     { in most cases, this function can be empty }
  351.     begin
  352.     end;
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361. {***********************************************************************************************}
  362.  
  363. { This is the main procedure through which all calls to the external module go.                    }
  364. { Main takes care of calling the right procedure with the right parameters depending on        }
  365. { the value of "selector".                                                                            }
  366. { You don't need to touch this procedure                                                            }
  367.     procedure mainMain (selector: integer; pb: ExtModulesParamBlockPtr);
  368.     begin
  369.         Startup(pb);
  370.         case selector of
  371.             kSetup: 
  372.                 begin
  373.                     pb^.requiredGlobals := 0;
  374.                     pb^.versionNumber := VERSIONNUMBER;
  375.                     if sizeof(extended) = 10 then
  376.                         pb^.codeType := CPU68noFPU
  377.                     else if sizeof(extended) = 12 then
  378.                         pb^.codeType := CPU68FPU
  379.                     else
  380.                         pb^.codeType := CPUPowerPC;
  381.  
  382.                     SetUp(pb^.moduleKind, pb^.name, pb^.requiredGlobals, pb);
  383.                 end;
  384.             progInitialize: 
  385.                 InitializeProg(pb);
  386.             progRun: 
  387.                 Run(pb);
  388.             kCleanUp: 
  389.                 CleanUp(pb);
  390.             otherwise
  391.         end;
  392.     end;
  393.  
  394.  
  395.  
  396. end.